sn-ENS neuron integration

for all current Baf53cre Ahr CTL/CKO data

#### sample info
## processed obj
##
#  CR7d    CTL  x4   CKO  x3 (x1 CKO failed)
#     20231220
#
#  Steady state   CTL  x3   CKO  x3
#     20240318
#
#
##
source("/Shared_win/projects/RNA_normal/analysis.10x.r")

load processed obj

merged.list <- list(CR7d.CTL_CKO=readRDS("I:/Shared_win/projects/20231220_10x_SZJ/analysis_plus/GEX0131.seur.Anno2.rds"),
                    Stst.CTL_CKO=readRDS("I:/Shared_win/projects/20240318_10x_SZJ/analysis_plus/GEX0430.seur.pure_Anno.rds")
                    )


merged.list
## $CR7d.CTL_CKO
## An object of class Seurat 
## 24081 features across 16520 samples within 1 assay 
## Active assay: RNA (24081 features, 1046 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
## 
## $Stst.CTL_CKO
## An object of class Seurat 
## 22951 features across 11932 samples within 1 assay 
## Active assay: RNA (22951 features, 1027 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
lapply(merged.list, function(x){unique(x$FB.info)})
## $CR7d.CTL_CKO
## [1] CTL.3 CTL.4 CTL.1 CKO.4 CKO.3 CKO.1 CTL.2
## Levels: CTL.1 CTL.2 CTL.3 CTL.4 CKO.1 CKO.3 CKO.4
## 
## $Stst.CTL_CKO
## [1] CTL.3 CKO.3 CTL.1 CTL.2 CKO.1 CKO.2
## Levels: CTL.1 CTL.2 CTL.3 CKO.1 CKO.2 CKO.3
lapply(merged.list, function(x){unique(x$cnt)})
## $CR7d.CTL_CKO
## [1] CTL CKO
## Levels: CTL CKO
## 
## $Stst.CTL_CKO
## [1] CTL CKO
## Levels: CTL CKO
lapply(merged.list, function(x){unique(x$rep)})
## $CR7d.CTL_CKO
## [1] "rep.3" "rep.4" "rep.1" "rep.2"
## 
## $Stst.CTL_CKO
## [1] "rep.3" "rep.1" "rep.2"

modify FB.info

rebuild cnt and rep as:

CR7d.CTL
CR7d.CKO

Stst.CTL
Stst.CKO

CR7d.CTL_CKO

#
merged.list$CR7d.CTL_CKO$sample <- as.character(merged.list$CR7d.CTL_CKO$FB.info)

merged.list$CR7d.CTL_CKO$sample <- as.vector(unlist(sapply(merged.list$CR7d.CTL_CKO$sample, function(x){
  gsub("CTL.","CR7d.CTL",x)
})))
merged.list$CR7d.CTL_CKO$sample <- as.vector(unlist(sapply(merged.list$CR7d.CTL_CKO$sample, function(x){
  gsub("CKO.","CR7d.CKO",x)
})))

merged.list$CR7d.CTL_CKO$cnt <- gsub("1$|2$|3$|4$","",as.character(merged.list$CR7d.CTL_CKO$sample))

merged.list$CR7d.CTL_CKO$sample[merged.list$CR7d.CTL_CKO$sample %in% "CR7d.CKO3"] <- "CR7d.CKO2"
merged.list$CR7d.CTL_CKO$sample[merged.list$CR7d.CTL_CKO$sample %in% "CR7d.CKO4"] <- "CR7d.CKO3"

merged.list$CR7d.CTL_CKO$rep <- paste0("rep",
                                       gsub("CR7d.CTL|CR7d.CKO","",as.character(merged.list$CR7d.CTL_CKO$sample)))

merged.list$CR7d.CTL_CKO$tissue <- "Ileum"
head(merged.list$CR7d.CTL_CKO$FB.info,10)
## AAACCCAAGAAGTCTA-1 AAACCCAAGACATACA-1 AAACCCAAGACTCATC-1 AAACCCAAGAGGTTTA-1 
##              CTL.3              CTL.3              CTL.4              CTL.1 
## AAACCCAAGGCTCTAT-1 AAACCCAAGGTCTTTG-1 AAACCCACAGATAAAC-1 AAACCCAGTATGAGAT-1 
##              CKO.4              CKO.3              CTL.4              CKO.1 
## AAACCCAGTATGGAAT-1 AAACCCAGTATGGGAC-1 
##              CTL.4              CTL.2 
## Levels: CTL.1 CTL.2 CTL.3 CTL.4 CKO.1 CKO.3 CKO.4
head(merged.list$CR7d.CTL_CKO$sample,10)
## AAACCCAAGAAGTCTA-1 AAACCCAAGACATACA-1 AAACCCAAGACTCATC-1 AAACCCAAGAGGTTTA-1 
##        "CR7d.CTL3"        "CR7d.CTL3"        "CR7d.CTL4"        "CR7d.CTL1" 
## AAACCCAAGGCTCTAT-1 AAACCCAAGGTCTTTG-1 AAACCCACAGATAAAC-1 AAACCCAGTATGAGAT-1 
##        "CR7d.CKO3"        "CR7d.CKO2"        "CR7d.CTL4"        "CR7d.CKO1" 
## AAACCCAGTATGGAAT-1 AAACCCAGTATGGGAC-1 
##        "CR7d.CTL4"        "CR7d.CTL2"
head(merged.list$CR7d.CTL_CKO$cnt,10)
## AAACCCAAGAAGTCTA-1 AAACCCAAGACATACA-1 AAACCCAAGACTCATC-1 AAACCCAAGAGGTTTA-1 
##         "CR7d.CTL"         "CR7d.CTL"         "CR7d.CTL"         "CR7d.CTL" 
## AAACCCAAGGCTCTAT-1 AAACCCAAGGTCTTTG-1 AAACCCACAGATAAAC-1 AAACCCAGTATGAGAT-1 
##         "CR7d.CKO"         "CR7d.CKO"         "CR7d.CTL"         "CR7d.CKO" 
## AAACCCAGTATGGAAT-1 AAACCCAGTATGGGAC-1 
##         "CR7d.CTL"         "CR7d.CTL"
head(merged.list$CR7d.CTL_CKO$rep,10)
## AAACCCAAGAAGTCTA-1 AAACCCAAGACATACA-1 AAACCCAAGACTCATC-1 AAACCCAAGAGGTTTA-1 
##             "rep3"             "rep3"             "rep4"             "rep1" 
## AAACCCAAGGCTCTAT-1 AAACCCAAGGTCTTTG-1 AAACCCACAGATAAAC-1 AAACCCAGTATGAGAT-1 
##             "rep3"             "rep2"             "rep4"             "rep1" 
## AAACCCAGTATGGAAT-1 AAACCCAGTATGGGAC-1 
##             "rep4"             "rep2"
head(merged.list$CR7d.CTL_CKO$tissue,10)
## AAACCCAAGAAGTCTA-1 AAACCCAAGACATACA-1 AAACCCAAGACTCATC-1 AAACCCAAGAGGTTTA-1 
##            "Ileum"            "Ileum"            "Ileum"            "Ileum" 
## AAACCCAAGGCTCTAT-1 AAACCCAAGGTCTTTG-1 AAACCCACAGATAAAC-1 AAACCCAGTATGAGAT-1 
##            "Ileum"            "Ileum"            "Ileum"            "Ileum" 
## AAACCCAGTATGGAAT-1 AAACCCAGTATGGGAC-1 
##            "Ileum"            "Ileum"

Stst.CTL_CKO

#
merged.list$Stst.CTL_CKO$sample <- as.character(merged.list$Stst.CTL_CKO$FB.info)

merged.list$Stst.CTL_CKO$sample <- as.vector(unlist(sapply(merged.list$Stst.CTL_CKO$sample, function(x){
  gsub("CTL.","Stst.CTL",x)
})))
merged.list$Stst.CTL_CKO$sample <- as.vector(unlist(sapply(merged.list$Stst.CTL_CKO$sample, function(x){
  gsub("CKO.","Stst.CKO",x)
})))

merged.list$Stst.CTL_CKO$cnt <- gsub("1$|2$|3$|4$","",as.character(merged.list$Stst.CTL_CKO$sample))

merged.list$Stst.CTL_CKO$rep <- paste0("rep",
                                       gsub("Stst.CTL|Stst.CKO","",as.character(merged.list$Stst.CTL_CKO$sample)))

merged.list$Stst.CTL_CKO$tissue <- "Ileum"
head(merged.list$Stst.CTL_CKO$FB.info,10)
## AAACCCAAGAATACAC-1 AAACCCAAGCAATAGT-1 AAACCCAAGGTGAGCT-1 AAACCCAAGTGGACGT-1 
##              CTL.3              CKO.3              CTL.1              CKO.3 
## AAACCCACAACGAGGT-1 AAACCCACAAGAGTAT-1 AAACCCACAATCAAGA-1 AAACCCACAGAGTCAG-1 
##              CTL.1              CKO.3              CTL.2              CKO.3 
## AAACCCAGTATCGATC-1 AAACCCATCAAGTTGC-1 
##              CTL.2              CKO.1 
## Levels: CTL.1 CTL.2 CTL.3 CKO.1 CKO.2 CKO.3
head(merged.list$Stst.CTL_CKO$sample,10)
## AAACCCAAGAATACAC-1 AAACCCAAGCAATAGT-1 AAACCCAAGGTGAGCT-1 AAACCCAAGTGGACGT-1 
##        "Stst.CTL3"        "Stst.CKO3"        "Stst.CTL1"        "Stst.CKO3" 
## AAACCCACAACGAGGT-1 AAACCCACAAGAGTAT-1 AAACCCACAATCAAGA-1 AAACCCACAGAGTCAG-1 
##        "Stst.CTL1"        "Stst.CKO3"        "Stst.CTL2"        "Stst.CKO3" 
## AAACCCAGTATCGATC-1 AAACCCATCAAGTTGC-1 
##        "Stst.CTL2"        "Stst.CKO1"
head(merged.list$Stst.CTL_CKO$cnt,10)
## AAACCCAAGAATACAC-1 AAACCCAAGCAATAGT-1 AAACCCAAGGTGAGCT-1 AAACCCAAGTGGACGT-1 
##         "Stst.CTL"         "Stst.CKO"         "Stst.CTL"         "Stst.CKO" 
## AAACCCACAACGAGGT-1 AAACCCACAAGAGTAT-1 AAACCCACAATCAAGA-1 AAACCCACAGAGTCAG-1 
##         "Stst.CTL"         "Stst.CKO"         "Stst.CTL"         "Stst.CKO" 
## AAACCCAGTATCGATC-1 AAACCCATCAAGTTGC-1 
##         "Stst.CTL"         "Stst.CKO"
head(merged.list$Stst.CTL_CKO$rep,10)
## AAACCCAAGAATACAC-1 AAACCCAAGCAATAGT-1 AAACCCAAGGTGAGCT-1 AAACCCAAGTGGACGT-1 
##             "rep3"             "rep3"             "rep1"             "rep3" 
## AAACCCACAACGAGGT-1 AAACCCACAAGAGTAT-1 AAACCCACAATCAAGA-1 AAACCCACAGAGTCAG-1 
##             "rep1"             "rep3"             "rep2"             "rep3" 
## AAACCCAGTATCGATC-1 AAACCCATCAAGTTGC-1 
##             "rep2"             "rep1"
head(merged.list$Stst.CTL_CKO$tissue,10)
## AAACCCAAGAATACAC-1 AAACCCAAGCAATAGT-1 AAACCCAAGGTGAGCT-1 AAACCCAAGTGGACGT-1 
##            "Ileum"            "Ileum"            "Ileum"            "Ileum" 
## AAACCCACAACGAGGT-1 AAACCCACAAGAGTAT-1 AAACCCACAATCAAGA-1 AAACCCACAGAGTCAG-1 
##            "Ileum"            "Ileum"            "Ileum"            "Ileum" 
## AAACCCAGTATCGATC-1 AAACCCATCAAGTTGC-1 
##            "Ileum"            "Ileum"

merge

all.batch.merged <- merge(merged.list[[1]],
                          merged.list[2:length(merged.list)],
                          merge.data=T)
## Warning in CheckDuplicateCellNames(object.list = objects): Some cell names are
## duplicated across objects provided. Renaming to enforce unique cell names.
all.batch.merged
## An object of class Seurat 
## 24355 features across 28452 samples within 1 assay 
## Active assay: RNA (24355 features, 0 variable features)
rm(merged.list)
gc()
##             used  (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   7790095 416.1   13105874   700.0   13105874   700.0
## Vcells 110003636 839.3 1401836563 10695.2 1710833604 13052.7
unique(all.batch.merged$cnt)
## [1] "CR7d.CTL" "CR7d.CKO" "Stst.CTL" "Stst.CKO"
unique(all.batch.merged$rep)
## [1] "rep3" "rep4" "rep1" "rep2"
all.batch.merged$cnt <- factor(as.character(all.batch.merged$cnt),
                               levels = c("Stst.CTL","Stst.CKO",
                                          "CR7d.CTL","CR7d.CKO"))
head(all.batch.merged$cnt)
## AAACCCAAGAAGTCTA-1_1 AAACCCAAGACATACA-1_1 AAACCCAAGACTCATC-1_1 
##             CR7d.CTL             CR7d.CTL             CR7d.CTL 
## AAACCCAAGAGGTTTA-1_1 AAACCCAAGGCTCTAT-1_1 AAACCCAAGGTCTTTG-1_1 
##             CR7d.CTL             CR7d.CKO             CR7d.CKO 
## Levels: Stst.CTL Stst.CKO CR7d.CTL CR7d.CKO
# modify orig.ident
all.batch.merged$orig.ident <- as.character(all.batch.merged$cnt)

all.batch.merged$orig.ident[all.batch.merged$orig.ident %in% c("Stst.CTL","Stst.CKO")] <- "Stst.CTL_CKO"
all.batch.merged$orig.ident[all.batch.merged$orig.ident %in% c("CR7d.CTL","CR7d.CKO")] <- "CR7d.CTL_CKO"


all.batch.merged$orig.ident <- factor(all.batch.merged$orig.ident,
                                      levels = c("Stst.CTL_CKO","CR7d.CTL_CKO"))

head(all.batch.merged$orig.ident)
## AAACCCAAGAAGTCTA-1_1 AAACCCAAGACATACA-1_1 AAACCCAAGACTCATC-1_1 
##         CR7d.CTL_CKO         CR7d.CTL_CKO         CR7d.CTL_CKO 
## AAACCCAAGAGGTTTA-1_1 AAACCCAAGGCTCTAT-1_1 AAACCCAAGGTCTTTG-1_1 
##         CR7d.CTL_CKO         CR7d.CTL_CKO         CR7d.CTL_CKO 
## Levels: Stst.CTL_CKO CR7d.CTL_CKO
table(all.batch.merged@meta.data[,c("orig.ident","cnt")])
##               cnt
## orig.ident     Stst.CTL Stst.CKO CR7d.CTL CR7d.CKO
##   Stst.CTL_CKO     5909     6023        0        0
##   CR7d.CTL_CKO        0        0    10280     6240
table(all.batch.merged@meta.data[,c("orig.ident","sample")])
##               sample
## orig.ident     CR7d.CKO1 CR7d.CKO2 CR7d.CKO3 CR7d.CTL1 CR7d.CTL2 CR7d.CTL3
##   Stst.CTL_CKO         0         0         0         0         0         0
##   CR7d.CTL_CKO      2365      1643      2232      2640      2507      2348
##               sample
## orig.ident     CR7d.CTL4 Stst.CKO1 Stst.CKO2 Stst.CKO3 Stst.CTL1 Stst.CTL2
##   Stst.CTL_CKO         0      1998      2218      1807      2122      1714
##   CR7d.CTL_CKO      2785         0         0         0         0         0
##               sample
## orig.ident     Stst.CTL3
##   Stst.CTL_CKO      2073
##   CR7d.CTL_CKO         0

Integration

#Idents(all.batch.merged) <- "cnt"
Idents(all.batch.merged) <- "orig.ident"

head(all.batch.merged@active.ident)
## AAACCCAAGAAGTCTA-1_1 AAACCCAAGACATACA-1_1 AAACCCAAGACTCATC-1_1 
##         CR7d.CTL_CKO         CR7d.CTL_CKO         CR7d.CTL_CKO 
## AAACCCAAGAGGTTTA-1_1 AAACCCAAGGCTCTAT-1_1 AAACCCAAGGTCTTTG-1_1 
##         CR7d.CTL_CKO         CR7d.CTL_CKO         CR7d.CTL_CKO 
## Levels: Stst.CTL_CKO CR7d.CTL_CKO
gc()
##             used  (Mb) gc trigger   (Mb)   max used    (Mb)
## Ncells   7791081 416.1   13105874  700.0   13105874   700.0
## Vcells 110006196 839.3 1121469251 8556.2 1710833604 13052.7
all.batch.merged
## An object of class Seurat 
## 24355 features across 28452 samples within 1 assay 
## Active assay: RNA (24355 features, 0 variable features)

SCT

batchcor.seu.list <- list()

#for(ttt in levels(all.batch.merged$cnt)){
#  ttt.seu <- subset(all.batch.merged, subset=cnt == ttt)
#  Idents(ttt.seu) <- "rep"
for(ttt in levels(all.batch.merged$orig.ident)){
  ttt.seu <- subset(all.batch.merged, subset=orig.ident == ttt)
  Idents(ttt.seu) <- "sample"
  ttt.seu <- SCTransform(ttt.seu, variable.features.n = 3000, vst.flavor="v2",
                         return.only.var.genes = FALSE)
  batchcor.seu.list[[ttt]] <- ttt.seu
}
## vst.flavor='v2' set, setting model to use fixed slope and exclude poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Total Step 1 genes: 17656
## Total overdispersed genes: 14190
## Excluding 3466 genes from Step 1 because they are not overdispersed.
## Variance stabilizing transformation of count matrix of size 19662 by 11932
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
## Setting estimate of  384 genes to inf as theta_mm/theta_mle < 1e-3
## # of step1 poisson genes (variance < mean): 0
## # of low mean genes (mean < 0.001): 2128
## Total # of Step1 poisson genes (theta=Inf; variance < mean): 410
## Total # of poisson genes (theta=Inf; variance < mean): 5429
## Calling offset model for all 5429 poisson genes
## Found 399 outliers - those will be ignored in fitting/regularization step
## Ignoring theta inf genes
## Replacing fit params for 5429 poisson genes by theta=Inf
## Setting min_variance based on median UMI:  0.04
## Second step: Get residuals using fitted parameters for 19662 genes
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |======================================================================| 100%
## Computing corrected count matrix for 19662 genes
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |======================================================================| 100%
## Calculating gene attributes
## Wall clock passed: Time difference of 1.36219 mins
## Determine variable features
## Place corrected count matrix in counts slot
## Centering data matrix
## Set default assay to SCT
## vst.flavor='v2' set, setting model to use fixed slope and exclude poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Total Step 1 genes: 18059
## Total overdispersed genes: 15302
## Excluding 2757 genes from Step 1 because they are not overdispersed.
## Variance stabilizing transformation of count matrix of size 20759 by 16520
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
## Setting estimate of  435 genes to inf as theta_mm/theta_mle < 1e-3
## # of step1 poisson genes (variance < mean): 0
## # of low mean genes (mean < 0.001): 2892
## Total # of Step1 poisson genes (theta=Inf; variance < mean): 465
## Total # of poisson genes (theta=Inf; variance < mean): 5437
## Calling offset model for all 5437 poisson genes
## Found 612 outliers - those will be ignored in fitting/regularization step
## Ignoring theta inf genes
## Replacing fit params for 5437 poisson genes by theta=Inf
## Setting min_variance based on median UMI:  0.04
## Second step: Get residuals using fitted parameters for 20759 genes
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |======================================================================| 100%
## Computing corrected count matrix for 20759 genes
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |======================================================================| 100%
## Calculating gene attributes
## Wall clock passed: Time difference of 1.872609 mins
## Determine variable features
## Place corrected count matrix in counts slot
## Centering data matrix
## Set default assay to SCT
gc()
##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   8092367  432.2   13105874   700.0   13105874   700.0
## Vcells 877719746 6696.5 2811167344 21447.6 2808264459 21425.4
rm(ttt.seu)
rm(all.batch.merged)
gc()
##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   8092284  432.2   13105874   700.0   13105874   700.0
## Vcells 780824647 5957.3 2248933876 17158.1 2808264459 21425.4
feature.raw <- SelectIntegrationFeatures(object.list = batchcor.seu.list, nfeatures = 3000)
head(feature.raw,300)
##   [1] "Cntn5"         "Kcnip4"        "Mgat4c"        "Cntnap2"      
##   [5] "Zfp804a"       "Klhl1"         "Kcnb2"         "Gal"          
##   [9] "Nrxn3"         "Pcdh9"         "Gm30613"       "Adgrg6"       
##  [13] "Lingo2"        "Asic2"         "Robo2"         "Ntng1"        
##  [17] "Rbfox1"        "Nrg1"          "Sgcz"          "Cpne4"        
##  [21] "Cdh18"         "Cmah"          "Gm32647"       "Prkg1"        
##  [25] "Schip1"        "Gm29521"       "Tmeff2"        "Ebf1"         
##  [29] "Ntrk3"         "Gm15261"       "Kcnq5"         "Cdh9"         
##  [33] "Brinp3"        "Egfem1"        "March1"        "Ano2"         
##  [37] "4930509J09Rik" "Cdh8"          "Nmu"           "Csmd3"        
##  [41] "Gm21847"       "Galntl6"       "Gpc6"          "Pdzrn4"       
##  [45] "Grik1"         "Nos1"          "Astn2"         "Dgkb"         
##  [49] "Kcnd2"         "4930555F03Rik" "Dgkg"          "Car10"        
##  [53] "Vwc2l"         "Ptprt"         "Zfp804b"       "Grm5"         
##  [57] "A330015K06Rik" "Dcc"           "Sema5a"        "Spock3"       
##  [61] "Wdr17"         "Unc5d"         "Chsy3"         "Sst"          
##  [65] "Vip"           "Nrg3"          "Pbx3"          "Csmd1"        
##  [69] "Fgf14"         "Gpr149"        "Tshz2"         "Zfhx3"        
##  [73] "Bnc2"          "Hs6st3"        "Trps1"         "Tafa1"        
##  [77] "Pde4d"         "Gpc5"          "Kcnh7"         "Tenm2"        
##  [81] "Gm45321"       "Lrrtm4"        "Serpini1"      "Cdh6"         
##  [85] "Nxph2"         "Robo1"         "Plxna4"        "Il1rapl1"     
##  [89] "Plcl1"         "4930432L08Rik" "Epha5"         "Gm20754"      
##  [93] "Rasgef1b"      "Tafa2"         "Dpp10"         "Clstn2"       
##  [97] "Gm15398"       "Gabrg3"        "Tac1"          "5530401A14Rik"
## [101] "Fam155a"       "Dach1"         "4930587E11Rik" "Ccbe1"        
## [105] "Gm42756"       "Lrrc4c"        "Auts2"         "Cadm2"        
## [109] "Nxph1"         "Grm7"          "Trhde"         "Col25a1"      
## [113] "Piezo2"        "9530059O14Rik" "Dock10"        "Cpa6"         
## [117] "Alcam"         "B230312C02Rik" "mt-Co3"        "1700111E14Rik"
## [121] "Tnr"           "Chrna9"        "Kctd8"         "Pde1a"        
## [125] "Skap1"         "4930445B16Rik" "Gm20658"       "Gm28376"      
## [129] "Gm20713"       "Penk"          "Gm15680"       "Gm36041"      
## [133] "Gm11342"       "Gm15584"       "Cntn3"         "Thsd7b"       
## [137] "Lrp1b"         "Cadps2"        "Efr3a"         "Itgb6"        
## [141] "Gm16638"       "Cntn4"         "Pcdh7"         "Grin3a"       
## [145] "Plcxd3"        "mt-Atp6"       "Ctnna3"        "Kcnab1"       
## [149] "Pcdh10"        "E130114P18Rik" "Ndst4"         "Dmd"          
## [153] "Brinp2"        "Cdc14a"        "Fbxw15"        "Gm12128"      
## [157] "Gm12239"       "mt-Co1"        "Sorcs1"        "Gm48742"      
## [161] "Gm11339"       "Caln1"         "Fhit"          "Cysltr2"      
## [165] "Gm48877"       "Tox"           "Myl1"          "Etv1"         
## [169] "D030068K23Rik" "C79798"        "Gm16226"       "Gm2516"       
## [173] "Pcdh11x"       "Galnt13"       "Cacna2d3"      "Slc44a5"      
## [177] "Galnt18"       "Specc1"        "Tmem132cos"    "Gm20631"      
## [181] "Rgs6"          "Rspo2"         "Fbxw24"        "4930428E07Rik"
## [185] "St6galnac3"    "Unc5c"         "Pcdh15"        "Gria1"        
## [189] "Ptprd"         "Synpr"         "Gna14"         "4930422I22Rik"
## [193] "Chrna7"        "Tmem132c"      "Grid2"         "Gfra1"        
## [197] "4930590L20Rik" "Slc4a4"        "Alk"           "Rtl4"         
## [201] "Ghr"           "Gm20642"       "Gm34623"       "Pde4b"        
## [205] "Dock2"         "Kirrel3"       "Gm28153"       "Pantr1"       
## [209] "Rarb"          "Fgf13"         "C130073E24Rik" "Mdga2"        
## [213] "Kctd16"        "Kcnj3"         "Kif26b"        "Nkain2"       
## [217] "Luzp2"         "Tmtc2"         "mt-Co2"        "Epha8"        
## [221] "Dgki"          "Slc35f4"       "Galnt17"       "Pde10a"       
## [225] "Gm10848"       "Sybu"          "Lama2"         "Trpm3"        
## [229] "Stxbp6"        "Naaladl2"      "Kcnt2"         "Kitl"         
## [233] "Gm4876"        "Lsamp"         "Gm43948"       "Pard3bos1"    
## [237] "Zfp536"        "Gm48094"       "Tenm3"         "Nlgn1"        
## [241] "Meis1"         "Nrg3os"        "Timp3"         "Pld5"         
## [245] "Kcnq3"         "Npas3"         "Prune2"        "9630028H03Rik"
## [249] "Htr4"          "Dsc3"          "Syn3"          "Fat3"         
## [253] "Airn"          "Cntnap5a"      "1700024B18Rik" "Dlc1"         
## [257] "Negr1"         "Gm6994"        "Frmd4b"        "Popdc3"       
## [261] "A830018L16Rik" "Eya1"          "Cux2"          "Grid1"        
## [265] "Pigk"          "mt-Nd2"        "mt-Nd4"        "Tenm4"        
## [269] "Pde7b"         "Adgrl3"        "Rmst"          "Zfp521"       
## [273] "Arhgap6"       "Tmem163"       "8030451O07Rik" "Gm11418"      
## [277] "Colec12"       "Dlgap2"        "Vcan"          "Mid1"         
## [281] "Ryr3"          "Col18a1"       "Pifo"          "Thsd4"        
## [285] "Arhgap24"      "Bmpr1b"        "Tenm1"         "Hmcn1"        
## [289] "Chrm2"         "Gm13269"       "Nek1"          "Nell1"        
## [293] "Cacnb2"        "Calcb"         "St18"          "Ank2"         
## [297] "Gm48321"       "Gm49678"       "Gm28154"       "Nfia"
DIG <- grep("^mt-|^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^Jun|^Fos|^AY|^AA|^AC|^AI|^AW|^BC|^Gm|^Hist|Lars2|Rik$|-ps",
            feature.raw,value = T)

# gender-relaged genes
GRG <- c("Xist","Tsix","Uty","Ddx3y","Eif2s3y","Kdm5d")

CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
DIG
##   [1] "Gm30613"       "Gm32647"       "Gm29521"       "Gm15261"      
##   [5] "4930509J09Rik" "Gm21847"       "4930555F03Rik" "A330015K06Rik"
##   [9] "Gm45321"       "4930432L08Rik" "Gm20754"       "Gm15398"      
##  [13] "5530401A14Rik" "4930587E11Rik" "Gm42756"       "9530059O14Rik"
##  [17] "B230312C02Rik" "mt-Co3"        "1700111E14Rik" "4930445B16Rik"
##  [21] "Gm20658"       "Gm28376"       "Gm20713"       "Gm15680"      
##  [25] "Gm36041"       "Gm11342"       "Gm15584"       "Gm16638"      
##  [29] "mt-Atp6"       "E130114P18Rik" "Gm12128"       "Gm12239"      
##  [33] "mt-Co1"        "Gm48742"       "Gm11339"       "Gm48877"      
##  [37] "D030068K23Rik" "Gm16226"       "Gm2516"        "Gm20631"      
##  [41] "4930428E07Rik" "4930422I22Rik" "4930590L20Rik" "Gm20642"      
##  [45] "Gm34623"       "Gm28153"       "C130073E24Rik" "mt-Co2"       
##  [49] "Gm10848"       "Gm4876"        "Gm43948"       "Gm48094"      
##  [53] "9630028H03Rik" "1700024B18Rik" "Gm6994"        "A830018L16Rik"
##  [57] "mt-Nd2"        "mt-Nd4"        "8030451O07Rik" "Gm11418"      
##  [61] "Gm13269"       "Gm48321"       "Gm49678"       "Gm28154"      
##  [65] "Gm15577"       "Gm32509"       "Gm42836"       "Gm45486"      
##  [69] "Gm27153"       "Gm11417"       "4930473D10Rik" "Gm49003"      
##  [73] "Gm12296"       "Gm48715"       "4930486F22Rik" "Gm28905"      
##  [77] "Gm26995"       "Gm42439"       "Gm20616"       "Gm16168"      
##  [81] "4930517O19Rik" "Gm21798"       "Gm26936"       "Gm11099"      
##  [85] "Gm13629"       "4930419G24Rik" "Gm27016"       "4933432K03Rik"
##  [89] "Gm48749"       "Gm4117"        "Gm45352"       "Gm47271"      
##  [93] "Gm37459"       "Gm3848"        "Gm11802"       "Gm45323"      
##  [97] "4930447N08Rik" "D030055H07Rik" "Gm26873"       "Gm15990"      
## [101] "Gm45341"       "Gm11844"       "Gm1992"        "Gm13963"      
## [105] "Gm19744"       "Gm26827"       "Gm28410"       "Gm49127"      
## [109] "Gm32335"       "Gm20641"       "Gm49959"       "2610307P16Rik"
## [113] "Gm35188"       "Gm48747"       "Gm11817"       "Gm48727"      
## [117] "Gm15336"       "Gm44151"       "4921534H16Rik" "Gm48093"      
## [121] "Gm38604"       "Gm14211"       "Gm48765"       "Gm16541"      
## [125] "Gm16599"       "Gm44071"       "Gm46124"       "A730090N16Rik"
## [129] "Gm12002"       "Gm50368"       "Gm12216"       "A230004M16Rik"
## [133] "Gm37679"       "Gm45645"       "Gm31698"       "9330175M20Rik"
## [137] "Gm36975"       "Gm11713"       "9630013A20Rik" "Gm16193"      
## [141] "Gm16685"       "4833422C13Rik" "Gm42413"       "Gm28501"      
## [145] "Gm15738"       "Gm16833"       "E330018M18Rik" "Gm13919"      
## [149] "B430010I23Rik" "B230307C23Rik" "Gm26672"       "Gm5441"       
## [153] "Gm40841"       "Gm47283"       "4930444A19Rik" "4933413L06Rik"
## [157] "Gm48228"       "4933424G05Rik" "Gm26954"       "Gm16070"      
## [161] "Gm15283"       "Gm26854"       "Gm31763"       "4921511C10Rik"
## [165] "Gm5535"        "Gm43569"       "2410021H03Rik" "Gm36198"      
## [169] "Gm16178"       "Gm12709"       "Gm19522"       "Gm33696"      
## [173] "4930545L23Rik" "Gm28375"       "mt-Cytb"       "Gm15563"      
## [177] "Gm16189"       "4930567K12Rik" "4930547E14Rik" "Gm41609"      
## [181] "Gm9968"        "mt-Nd3"        "4930484I04Rik" "A230009B12Rik"
## [185] "4930570B17Rik" "Gm4258"        "Gm11149"       "Gm33228"      
## [189] "Gm11437"       "Gm49961"       "5730522E02Rik" "Gm10550"      
## [193] "Gm43376"       "Gm36146"       "Gm27151"       "Gm41505"      
## [197] "3110021N24Rik" "4930595D18Rik" "Gm15337"       "Gm42556"      
## [201] "Gm15614"       "Gm11536"       "4930511M06Rik" "Gm15886"      
## [205] "4930447A16Rik" "Gm47027"       "Gm15155"       "A330008L17Rik"
## [209] "Gm34030"       "Gm50270"       "Gm5089"        "1810062O18Rik"
## [213] "4930592A05Rik" "Gm28750"       "Gm15496"       "Gm48383"      
## [217] "Gm47071"       "Gm14798"       "Gm49226"       "Gm48342"      
## [221] "1810034E14Rik" "Gm16183"       "Gm26911"       "Gm2694"       
## [225] "Gm16759"       "mt-Nd1"        "Gm15904"       "Gm34299"      
## [229] "Rps29"         "Gm16093"       "Gm19605"       "A930002C04Rik"
## [233] "Gm48512"       "C330002G04Rik" "Gm39214"       "mt-Nd5"       
## [237] "4933408N05Rik" "Gm32014"       "Gm16036"       "A930037H05Rik"
## [241] "Gm49179"       "Gm39090"       "Gm37240"       "Gm37229"      
## [245] "Rpl37"         "Gm10563"       "Gm11457"       "Rps27a"       
## [249] "9630014M24Rik" "1700109K24Rik" "Gm28198"       "A330102I10Rik"
## [253] "Gm2087"        "Gm43000"       "2610204G07Rik" "2700054A10Rik"
## [257] "4930414N06Rik" "4932443L11Rik" "4930532I03Rik" "Rpl39"        
## [261] "Gm17231"       "Gm19710"       "4930512B01Rik" "Gm43915"      
## [265] "Gm31683"       "Gm29530"       "Gm17178"       "Gm20470"      
## [269] "2810432F15Rik" "Traf3ip2"      "Gm26918"       "Gm17545"      
## [273] "Gm42418"       "Gm10135"       "Gm38251"       "Gm26542"      
## [277] "Gm26691"       "1700037C18Rik" "Gm38304"       "C630031E19Rik"
## [281] "A830019P07Rik" "Gm14858"       "4930527B05Rik" "Gm7580"       
## [285] "Gm42477"       "9330188P03Rik" "Gm12353"       "Gm27241"      
## [289] "Gm19461"       "6530409C15Rik" "Gm13402"       "Gm16958"      
## [293] "Gm50013"       "Gm37640"       "Gm12592"       "Gm50163"      
## [297] "Gm10785"       "1700085D07Rik" "AC154200.1"    "Gm38134"      
## [301] "Gm50152"       "4831440D22Rik" "Gm46367"       "A330032B11Rik"
## [305] "Gm17058"       "Gm16066"       "Gm15478"       "Gm15700"      
## [309] "Gm43713"       "Gm49417"       "1700084C06Rik" "Gm20149"      
## [313] "Rps6ka3"       "Gm50232"       "Gm17227"       "Gm42609"      
## [317] "D830024N08Rik" "Gm30054"       "4930557K07Rik" "Gm26759"      
## [321] "Gm17036"       "1600010M07Rik" "Rpl24"         "Gm20536"      
## [325] "Gm37571"       "Rps21"         "Gm48677"       "Gm15489"      
## [329] "Gm49521"       "3110082I17Rik" "Gm5820"        "Gm10544"      
## [333] "Gm19667"       "Gm47153"       "Gm26724"       "Gm26658"      
## [337] "Gm39326"       "Gm10629"       "Gm10658"       "AC163638.2"   
## [341] "2610027K06Rik" "A230006K03Rik" "Rpl26"         "Rpl37a"       
## [345] "Gm47167"       "Gm48708"       "Gm50020"       "1110046J04Rik"
## [349] "Gm15396"       "Rpsa"          "Gm16043"       "Rps28"        
## [353] "Gm10130"       "Gm43192"       "A230072E10Rik" "2610020C07Rik"
## [357] "D830036C21Rik" "4933400F21Rik" "Gm37768"       "2610316D01Rik"
## [361] "Dnajc15"       "B020031H02Rik" "Gm10791"       "Rpl38"        
## [365] "Rplp1"         "6430590A07Rik" "4933433G19Rik" "mt-Nd6"       
## [369] "4930583P06Rik" "Rps23"         "4930455C13Rik" "Gm44686"      
## [373] "Rplp2"         "Hsp90ab1"      "Dnajc19"       "5033428I22Rik"
## [377] "Gm31323"       "4930556N13Rik" "2810004N23Rik" "Gm29114"      
## [381] "Rps24"         "Traf3"         "Rpl36"         "Rpl27"        
## [385] "Gm19466"       "Gm47134"       "Gm20404"       "Gm17529"      
## [389] "Gm10353"       "1700126G02Rik" "2700069I18Rik" "BC052040"     
## [393] "Gm15787"       "Gm16794"       "C130071C03Rik" "Gm15270"      
## [397] "Gm26979"       "Gm50250"       "Gm16552"       "A630089N07Rik"
## [401] "Gm48855"       "Gm15952"       "Gm41764"       "Gm46440"      
## [405] "1700110K17Rik" "Gm34868"       "A530058N18Rik" "4930506C21Rik"
## [409] "Gm28085"       "Rps20"         "Rpl30"         "Gm15243"      
## [413] "Rps15a"        "4921511I17Rik" "E530011L22Rik" "Gm11659"      
## [417] "Gm15893"       "Gm49915"       "1700012D14Rik" "2210408I21Rik"
## [421] "Gm14372"       "5330417C22Rik" "Rps12"         "4930480K15Rik"
## [425] "Gm47350"       "4930429F24Rik" "Rpl23"         "Gm50138"      
## [429] "0610043K17Rik" "Gm11906"       "Gm26788"       "Gm48239"      
## [433] "AW554918"      "B230216N24Rik" "Rpl32"         "6530437J22Rik"
## [437] "Gm16279"       "Gm28379"       "Rpl35a"        "Gm49189"      
## [441] "Gm45509"       "Gm3294"        "Gm45435"       "Gm47551"      
## [445] "Gm13187"       "9330182L06Rik" "4930565N06Rik" "Gm17494"      
## [449] "1700087M22Rik" "Gm17477"       "Gm26749"       "Rpl27a"       
## [453] "Rps9"          "4930473A02Rik" "4930556J24Rik" "Gm16185"      
## [457] "Rps14"         "4930586N03Rik" "4930549G23Rik" "Rpl13a"       
## [461] "Rps26"         "Gm26981"       "Rps27"         "C530008M17Rik"
## [465] "Rps3a1"        "Gm1604a"       "Rpl8"          "Rpl11"        
## [469] "Gm32743"       "Gm37982"       "Rpl10a"        "Rpl3"         
## [473] "Rps6ka5"       "E130307A14Rik" "AW551984"      "Rps8"         
## [477] "A230056P14Rik" "Gm11250"       "6330403K07Rik" "Gm26531"      
## [481] "Gm26673"       "Rpl6"          "Gm12473"       "E230016M11Rik"
## [485] "Dnajc6"        "Gm44649"       "1700080N15Rik" "6430548M08Rik"
## [489] "Rpl5"          "Gm13748"       "Gmcl1"         "4933421D24Rik"
## [493] "Rpl12"         "Trafd1"        "Gm47460"       "Gm20721"      
## [497] "2610035D17Rik" "Gm36992"       "Gm41804"       "Gm10125"      
## [501] "Rpl19"         "Rps17"         "Gm19938"       "Dnaja3"       
## [505] "Gm48678"       "Gm16083"       "Gm49336"       "Gmeb2"        
## [509] "1500004A13Rik" "Rpl21"         "D430019H16Rik" "Gm12940"      
## [513] "B130055M24Rik" "Gm47802"       "4933406I18Rik" "2700081O15Rik"
## [517] "Gm47827"       "Gm30648"       "Gm15632"       "D430042O09Rik"
## [521] "Gm28370"       "Gm32250"       "Gm16549"       "1110002L01Rik"
## [525] "Rps4x"         "Dnajc24"       "Rpl41"         "2510009E07Rik"
## [529] "Rpl13"         "1110059E24Rik" "Gm43388"       "9530026P05Rik"
## [533] "1110019D14Rik" "mt-Atp8"       "mt-Nd4l"       "4931403E22Rik"
## [537] "Gm49191"       "Gm49171"       "Gm45600"       "Gm13708"      
## [541] "1700094D03Rik" "AC149090.1"    "Trank1"        "Rpl22l1"      
## [545] "Gm26879"       "Rpl34"         "4930578M01Rik" "Rpl35"        
## [549] "AC150683.1"    "Rps18"         "Gm12227"       "Gm16351"      
## [553] "Ndufb1-ps"     "Rps3"          "Gm44129"       "Hsp90aa1"     
## [557] "5430427M07Rik" "Rpl18"         "Gm15764"       "9330185C12Rik"
## [561] "Rps10"         "Dnajc7"        "6720427I07Rik" "Dnajc3"       
## [565] "AW549877"      "Rps19"         "Gm29237"       "Gm47889"      
## [569] "Gm17202"       "Gm27008"       "Rps7"          "Gm10076"      
## [573] "Gm12298"       "Gm42788"       "Gm47661"       "Rps13"        
## [577] "4930599N23Rik" "Gm12224"       "1600020E01Rik" "D630045J12Rik"
## [581] "Gm26618"       "Gm20456"       "Gm3227"        "AI480526"     
## [585] "Gm30094"       "Rpl7a"         "Gm11867"       "Gm10974"      
## [589] "Gm11465"       "2610037D02Rik" "Gm11714"       "Rpl18a"       
## [593] "Gm16068"       "Gm36279"       "BC048559"      "2010300C02Rik"
## [597] "2210408F21Rik" "5730480H06Rik" "3110056K07Rik" "Rpl9"         
## [601] "2700049A03Rik" "Hspa5"         "Hist3h2ba"     "Gm45502"      
## [605] "Gm36363"       "Rps6"          "1810026B05Rik" "Gm28294"
CC_gene
##  [1] "Mcm5"     "Pcna"     "Tyms"     "Fen1"     "Mcm7"     "Mcm4"    
##  [7] "Rrm1"     "Ung"      "Gins2"    "Mcm6"     "Cdca7"    "Dtl"     
## [13] "Prim1"    "Uhrf1"    "Cenpu"    "Hells"    "Rfc2"     "Polr1b"  
## [19] "Nasp"     "Rad51ap1" "Gmnn"     "Wdr76"    "Slbp"     "Ccne2"   
## [25] "Ubr7"     "Pold3"    "Msh2"     "Atad2"    "Rad51"    "Rrm2"    
## [31] "Cdc45"    "Cdc6"     "Exo1"     "Tipin"    "Dscc1"    "Blm"     
## [37] "Casp8ap2" "Usp1"     "Clspn"    "Pola1"    "Chaf1b"   "Mrpl36"  
## [43] "E2f8"     "Hmgb2"    "Cdk1"     "Nusap1"   "Ube2c"    "Birc5"   
## [49] "Tpx2"     "Top2a"    "Ndc80"    "Cks2"     "Nuf2"     "Cks1b"   
## [55] "Mki67"    "Tmpo"     "Cenpf"    "Tacc3"    "Pimreg"   "Smc4"    
## [61] "Ccnb2"    "Ckap2l"   "Ckap2"    "Aurkb"    "Bub1"     "Kif11"   
## [67] "Anp32e"   "Tubb4b"   "Gtse1"    "Kif20b"   "Hjurp"    "Cdca3"   
## [73] "Jpt1"     "Cdc20"    "Ttk"      "Cdc25c"   "Kif2c"    "Rangap1" 
## [79] "Ncapd2"   "Dlgap5"   "Cdca2"    "Cdca8"    "Ect2"     "Kif23"   
## [85] "Hmmr"     "Aurka"    "Psrc1"    "Anln"     "Lbr"      "Ckap5"   
## [91] "Cenpe"    "Ctcf"     "Nek2"     "G2e3"     "Gas2l3"   "Cbx5"    
## [97] "Cenpa"
features.filt <- setdiff(feature.raw, c(DIG,CC_gene,GRG))
length(features.filt)
## [1] 2386
head(features.filt,300)
##   [1] "Cntn5"      "Kcnip4"     "Mgat4c"     "Cntnap2"    "Zfp804a"   
##   [6] "Klhl1"      "Kcnb2"      "Gal"        "Nrxn3"      "Pcdh9"     
##  [11] "Adgrg6"     "Lingo2"     "Asic2"      "Robo2"      "Ntng1"     
##  [16] "Rbfox1"     "Nrg1"       "Sgcz"       "Cpne4"      "Cdh18"     
##  [21] "Cmah"       "Prkg1"      "Schip1"     "Tmeff2"     "Ebf1"      
##  [26] "Ntrk3"      "Kcnq5"      "Cdh9"       "Brinp3"     "Egfem1"    
##  [31] "March1"     "Ano2"       "Cdh8"       "Nmu"        "Csmd3"     
##  [36] "Galntl6"    "Gpc6"       "Pdzrn4"     "Grik1"      "Nos1"      
##  [41] "Astn2"      "Dgkb"       "Kcnd2"      "Dgkg"       "Car10"     
##  [46] "Vwc2l"      "Ptprt"      "Zfp804b"    "Grm5"       "Dcc"       
##  [51] "Sema5a"     "Spock3"     "Wdr17"      "Unc5d"      "Chsy3"     
##  [56] "Sst"        "Vip"        "Nrg3"       "Pbx3"       "Csmd1"     
##  [61] "Fgf14"      "Gpr149"     "Tshz2"      "Zfhx3"      "Bnc2"      
##  [66] "Hs6st3"     "Trps1"      "Tafa1"      "Pde4d"      "Gpc5"      
##  [71] "Kcnh7"      "Tenm2"      "Lrrtm4"     "Serpini1"   "Cdh6"      
##  [76] "Nxph2"      "Robo1"      "Plxna4"     "Il1rapl1"   "Plcl1"     
##  [81] "Epha5"      "Rasgef1b"   "Tafa2"      "Dpp10"      "Clstn2"    
##  [86] "Gabrg3"     "Tac1"       "Fam155a"    "Dach1"      "Ccbe1"     
##  [91] "Lrrc4c"     "Auts2"      "Cadm2"      "Nxph1"      "Grm7"      
##  [96] "Trhde"      "Col25a1"    "Piezo2"     "Dock10"     "Cpa6"      
## [101] "Alcam"      "Tnr"        "Chrna9"     "Kctd8"      "Pde1a"     
## [106] "Skap1"      "Penk"       "Cntn3"      "Thsd7b"     "Lrp1b"     
## [111] "Cadps2"     "Efr3a"      "Itgb6"      "Cntn4"      "Pcdh7"     
## [116] "Grin3a"     "Plcxd3"     "Ctnna3"     "Kcnab1"     "Pcdh10"    
## [121] "Ndst4"      "Dmd"        "Brinp2"     "Cdc14a"     "Fbxw15"    
## [126] "Sorcs1"     "Caln1"      "Fhit"       "Cysltr2"    "Tox"       
## [131] "Myl1"       "Etv1"       "C79798"     "Pcdh11x"    "Galnt13"   
## [136] "Cacna2d3"   "Slc44a5"    "Galnt18"    "Specc1"     "Tmem132cos"
## [141] "Rgs6"       "Rspo2"      "Fbxw24"     "St6galnac3" "Unc5c"     
## [146] "Pcdh15"     "Gria1"      "Ptprd"      "Synpr"      "Gna14"     
## [151] "Chrna7"     "Tmem132c"   "Grid2"      "Gfra1"      "Slc4a4"    
## [156] "Alk"        "Rtl4"       "Ghr"        "Pde4b"      "Dock2"     
## [161] "Kirrel3"    "Pantr1"     "Rarb"       "Fgf13"      "Mdga2"     
## [166] "Kctd16"     "Kcnj3"      "Kif26b"     "Nkain2"     "Luzp2"     
## [171] "Tmtc2"      "Epha8"      "Dgki"       "Slc35f4"    "Galnt17"   
## [176] "Pde10a"     "Sybu"       "Lama2"      "Trpm3"      "Stxbp6"    
## [181] "Naaladl2"   "Kcnt2"      "Kitl"       "Lsamp"      "Pard3bos1" 
## [186] "Zfp536"     "Tenm3"      "Nlgn1"      "Meis1"      "Nrg3os"    
## [191] "Timp3"      "Pld5"       "Kcnq3"      "Npas3"      "Prune2"    
## [196] "Htr4"       "Dsc3"       "Syn3"       "Fat3"       "Airn"      
## [201] "Cntnap5a"   "Dlc1"       "Negr1"      "Frmd4b"     "Popdc3"    
## [206] "Eya1"       "Cux2"       "Grid1"      "Pigk"       "Tenm4"     
## [211] "Pde7b"      "Adgrl3"     "Rmst"       "Zfp521"     "Arhgap6"   
## [216] "Tmem163"    "Colec12"    "Dlgap2"     "Vcan"       "Mid1"      
## [221] "Ryr3"       "Col18a1"    "Pifo"       "Thsd4"      "Arhgap24"  
## [226] "Bmpr1b"     "Tenm1"      "Hmcn1"      "Chrm2"      "Nek1"      
## [231] "Nell1"      "Cacnb2"     "Calcb"      "St18"       "Ank2"      
## [236] "Nfia"       "Sgcd"       "Synpo2"     "Khdrbs2"    "Dach2"     
## [241] "Sdk1"       "Nell1os"    "Htr2b"      "Adcy2"      "Syt6"      
## [246] "Olfm3"      "Mef2c"      "Slc26a4"    "Creb5"      "Plcb1"     
## [251] "Lrrc7"      "Fbxl7"      "Cdh20"      "Ptprz1"     "Cntnap5b"  
## [256] "Arpp21"     "Scg2"       "Man1a"      "Xylt1"      "Adamtsl1"  
## [261] "Zbbx"       "Cdh11"      "Edil3"      "Mir99ahg"   "Ammecr1"   
## [266] "Mast4"      "Arhgap15"   "Rims1"      "Eda"        "Oprk1"     
## [271] "Inpp4b"     "Tcf7l2"     "Hs3st2"     "Nfib"       "Ngfr"      
## [276] "Flrt2"      "Cpne8"      "Samd4"      "Kcnk2"      "Grm3"      
## [281] "Unc13c"     "Cdh10"      "Gda"        "Ano5"       "Plpp4"     
## [286] "Tafa5"      "Col5a2"     "Gabrb1"     "Prkca"      "Frmpd4"    
## [291] "Nxn"        "Gucy1a2"    "Usp6nl"     "Itgbl1"     "Bcl2"      
## [296] "Hccs"       "Igf2r"      "Thsd7a"     "Casz1"      "Ntrk2"
length(feature.raw)
## [1] 3000
length(features.filt)
## [1] 2386
GRG %in% feature.raw
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
GRG %in% features.filt
## [1] FALSE FALSE FALSE FALSE FALSE FALSE

anchors

all.anchors <- FindIntegrationAnchors(object.list = batchcor.seu.list,
                                      dims = 1:50,
                                      anchor.features = features.filt)
## Scaling features for provided objects
## Finding all pairwise anchors
## Running CCA
## Merging objects
## Finding neighborhoods
## Finding anchors
##  Found 37773 anchors
## Filtering anchors
##  Retained 9109 anchors
all.anchors
## An AnchorSet object containing 18218 anchors between 2 Seurat objects 
##  This can be used as input to IntegrateData.

Integration

all.anchors <- IntegrateData(anchorset = all.anchors, dims = 1:50)
## Merging dataset 1 into 2
## Extracting anchors for merged samples
## Finding integration vectors
## Finding integration vector weights
## Integrating data
all.anchors
## An object of class Seurat 
## 47749 features across 28452 samples within 3 assays 
## Active assay: integrated (2386 features, 2386 variable features)
##  2 other assays present: RNA, SCT
all.anchors@assays$SCT@SCTModel.list
## $model1
## An sctransform model.
##   Model formula:  y ~ log_umi 
##   Parameters stored for 19662 features, 11932 cells.
## 
## $model1.1
## An sctransform model.
##   Model formula:  y ~ log_umi 
##   Parameters stored for 20759 features, 16520 cells.
#rm(batchcor.seu.list)
gc()
##              used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells    8167903  436.3   13105874   700.0   13105874   700.0
## Vcells 1089076968 8309.0 2248933876 17158.1 2808264459 21425.4
#head(all.anchors$orig.ident)
#unique(all.anchors$orig.ident)
unique(all.anchors$cnt)
## [1] "Stst.CTL" "Stst.CKO" "CR7d.CTL" "CR7d.CKO"
#
all.anchors$orig.ident <- factor(all.anchors$orig.ident,
                                 levels = c("Stst.CTL_CKO","CR7d.CTL_CKO"))

# 
#table(all.anchors@meta.data[,c("cnt","condition")])
table(all.anchors@meta.data[,c("cnt","orig.ident")])
##           orig.ident
## cnt        Stst.CTL_CKO CR7d.CTL_CKO
##   CR7d.CKO            0         6240
##   CR7d.CTL            0        10280
##   Stst.CKO         6023            0
##   Stst.CTL         5909            0
VlnPlot(all.anchors, features = GRG, group.by = "sample", assay = "RNA")

all males

clustering

all.anchors <- ScaleData(object = all.anchors, verbose = TRUE, 
                         vars.to.regress = c("percent.mt","percent.rb","nCount_RNA"))
## Regressing out percent.mt, percent.rb, nCount_RNA
## Centering and scaling data matrix
# sample condition to regress
#   cnt ?
#   orig.idnet ?
#   condition ?
#   sex ?
gc()
##              used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells    8290322  442.8   13105874   700.0   13105874   700.0
## Vcells 1157659183 8832.3 2248933876 17158.1 2808264459 21425.4
length(all.anchors@assays$integrated@var.features)
## [1] 2386
all.anchors <- RunPCA(all.anchors, do.print = TRUE,
                      features = all.anchors@assays$integrated@var.features,
                      seed.use = 133,
                      npcs = 100,
                      #ndims.print = 1,
                      verbose = T)
## PC_ 1 
## Positive:  Nrg3, Grid2, Lrrtm4, Pcdh15, Trpm3, Ndst4, Cadm1, Kcnq3, Frmpd4, Epha5 
##     Cacnb2, Slit2, Kcnc2, Nrg3os, Kcnd2, Nkain2, Fstl5, Syn3, Magi1, Galntl6 
##     Ctnna3, Nav3, Lrrc4c, Kcnab1, Plce1, Kcnma1, Nos1, Pde4b, Lrp1b, Dmd 
## Negative:  Ntrk3, Ano2, Robo2, Cdh8, Nrxn3, Tmeff2, Cpne4, Clstn2, Myl1, Plxna4 
##     Pdzrn4, Adgrg6, Mgat4c, Zfp804a, Gpr149, Pcdh10, Ccbe1, Cysltr2, Cdh6, Dgkg 
##     Cux2, Spock3, Astn2, Grin3a, Iqgap2, Cntn5, Itgb6, Arhgap6, Kcnb2, Pcdh9 
## PC_ 2 
## Positive:  Bnc2, Rbfox1, Ptprt, Gpc6, Tshz2, Grik1, Tafa1, Mdga2, Tox, St6galnac3 
##     Adgrb3, Frmd4b, Pcdh7, Brinp2, Fbxw15, Cdc14a, Plcxd3, Agbl4, Negr1, Pde4b 
##     Pld5, Oprk1, Fbxw24, Dlgap2, Specc1, Chat, Xylt1, Unc5c, Dock2, Adamtsl1 
## Negative:  Auts2, Nos1, Dgkb, Alcam, Etv1, Fam155a, Cadps2, Gfra1, Kcnq5, Egfem1 
##     Kcnt2, Asic2, Kcnab1, Dach1, Schip1, Plekha5, Rgs6, Ank2, Cmah, Epha5 
##     Cntnap5a, Nav3, Stxbp6, Hs6st3, Tmem108, Creb5, Ncam2, Il1rapl1, Ablim2, Ebf1 
## PC_ 3 
## Positive:  Kcnip4, Cdh18, Csmd3, Klhl1, Kctd8, Pbx3, Cadm2, Gabrg3, Htr4, Cntn3 
##     Dlc1, Pde4d, Meis1, March1, Serpini1, C79798, Sema5a, Skap1, Car10, Khdrbs2 
##     Prkg1, Edil3, Zfhx3, Vwc2l, Plcl1, Gpc6, Dmd, Cdh9, Nrp2, Pakap 
## Negative:  Sgcd, Ptprg, Adgrg6, Nfia, Fgf13, Filip1, Nfib, Cysltr2, Sema6d, Gpr149 
##     Ccbe1, Malat1, Slc4a4, Nos1, Nmu, Cntnap5a, Grin3a, Ano2, Dapk2, Tshz2 
##     Airn, Dgkg, Zfp804a, Cbln2, Itgb6, Efr3a, Cpne4, Zfp521, Ryr2, Zfp536 
## PC_ 4 
## Positive:  Lingo2, Dock10, Prkg1, Kcnt2, Ndst4, Fgf13, Epha5, Kcnip4, Ctnna3, Gda 
##     Dmd, Lrrtm4, Tac1, Cntn5, Lrrc4c, Nxph1, Sorcs1, Chl1, Kctd8, Thsd4 
##     Ptprz1, Grem2, Ank2, Lama2, Hs3st2, Rgs6, Hgf, Rora, Lrrc7, Galntl6 
## Negative:  Chsy3, Trhde, March1, Ebf1, Trps1, Cntn4, Gal, Enox1, Nrg1, Trpm3 
##     Col18a1, Ntng1, Zmat4, Sdk1, Cpa6, Dcc, Ccser1, Tenm4, Shisa6, Plcxd3 
##     Csmd1, Npas3, Nkain3, Kcnh7, Sctr, Prune2, Tenm1, Sez6l, Sema5a, Galnt17 
## PC_ 5 
## Positive:  Nrg1, Trhde, Ptprd, Cntn4, Lsamp, Egfem1, Rmst, Cntn3, Adgrl2, Kcnd2 
##     Sgcz, Csmd3, Gal, Cntn5, Luzp2, Trps1, Ebf1, Nav2, Cpa6, Zmat4 
##     Hs6st3, Csmd1, Asic2, Astn2, Moxd1, Col18a1, Nkain3, Lrp1b, Rasgrf2, Kcnip4 
## Negative:  Dgkb, Klhl1, Vwc2l, Rasgef1b, Pbx3, Il1rapl1, Alk, Cdh9, Mgat4c, Alcam 
##     Zfhx3, Dpp10, Vcan, Sema5a, Fam155a, Galnt18, Auts2, Galnt13, Zbbx, Olfr78 
##     Thsd7b, C79798, Pcdh7, Scgn, Nek1, Serpini1, P3h2, Lncbate1, Stard13, Gucy1a2
DimHeatmap(all.anchors, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

ElbowPlot(all.anchors, ndims = 100)

ElbowPlot(all.anchors, ndims = 50)

Integrated

all.anchors@meta.data[,grep("snn|pANN|cluster|sort|preAnno|snn|Doublet",colnames(all.anchors@meta.data))] <- NULL

head(all.anchors@meta.data)
##                        orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAATACAC-1_2 Stst.CTL_CKO       2282         1290 0.13146363  0.1314636
## AAACCCAAGCAATAGT-1_2 Stst.CTL_CKO       1535         1011 0.13029316  0.5863192
## AAACCCAAGGTGAGCT-1_2 Stst.CTL_CKO       1049          737 0.00000000  0.4766444
## AAACCCAAGTGGACGT-1_2 Stst.CTL_CKO       4119         1907 0.09711095  0.4855547
## AAACCCACAACGAGGT-1_2 Stst.CTL_CKO       1591         1046 0.69138906  0.4399749
## AAACCCACAAGAGTAT-1_2 Stst.CTL_CKO       2106         1250 0.28490028  0.2849003
##                      FB.info     S.Score    G2M.Score Phase      cnt Anno1
## AAACCCAAGAATACAC-1_2   CTL.3 -0.00110357  0.012058898   G2M Stst.CTL  EMN1
## AAACCCAAGCAATAGT-1_2   CKO.3  0.02907907 -0.015558699     S Stst.CKO  IMN1
## AAACCCAAGGTGAGCT-1_2   CTL.1 -0.01100811 -0.008015087    G1 Stst.CTL  EMN1
## AAACCCAAGTGGACGT-1_2   CKO.3  0.01407052 -0.028388278     S Stst.CKO IPAN1
## AAACCCACAACGAGGT-1_2   CTL.1 -0.02143685  0.005086498   G2M Stst.CTL  EMN1
## AAACCCACAAGAGTAT-1_2   CKO.3  0.02965845 -0.009057774     S Stst.CKO  IMN3
##                        Anno2  rep Anno3    sample tissue nCount_SCT
## AAACCCAAGAATACAC-1_2    EMN1 rep3  <NA> Stst.CTL3  Ileum       1835
## AAACCCAAGCAATAGT-1_2    IMN1 rep3  <NA> Stst.CKO3  Ileum       1504
## AAACCCAAGGTGAGCT-1_2    EMN1 rep1  <NA> Stst.CTL1  Ileum       1159
## AAACCCAAGTGGACGT-1_2 IPAN1.1 rep3  <NA> Stst.CKO3  Ileum       1898
## AAACCCACAACGAGGT-1_2    EMN1 rep1  <NA> Stst.CTL1  Ileum       1550
## AAACCCACAAGAGTAT-1_2    IMN3 rep3  <NA> Stst.CKO3  Ileum       1799
##                      nFeature_SCT
## AAACCCAAGAATACAC-1_2         1263
## AAACCCAAGCAATAGT-1_2          989
## AAACCCAAGGTGAGCT-1_2          726
## AAACCCAAGTGGACGT-1_2         1223
## AAACCCACAACGAGGT-1_2         1027
## AAACCCACAAGAGTAT-1_2         1219
DefaultAssay(all.anchors) <- "integrated"
PCsct <- 1:21
all.anchors <- FindNeighbors(all.anchors, k.param = 20, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)
## Computing nearest neighbor graph
## Computing SNN
all.anchors <- FindClusters(all.anchors, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution =2, reduction = 'pca', verbose = T)
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Warning: The following arguments are not used: dims.use, save.SNN, reduction
## Suggested parameter: dims instead of dims.use
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 28452
## Number of edges: 1093172
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8385
## Number of communities: 32
## Elapsed time: 5 seconds
all.anchors <- RunTSNE(object = all.anchors, assay = "integrated", seed.use = 233, dims = PCsct, complexity=100)
all.anchors <- RunUMAP(object = all.anchors, assay = "integrated", seed.use = 888, dims = PCsct, n.neighbors = 20, min.dist = 0.3)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 16:27:40 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:27:40 Read 28452 rows and found 21 numeric columns
## 16:27:40 Using Annoy for neighbor search, n_neighbors = 20
## 16:27:40 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:27:44 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpCcsu6f\file760018e64cd4
## 16:27:44 Searching Annoy index using 1 thread, search_k = 2000
## 16:27:50 Annoy recall = 100%
## 16:27:51 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 16:27:53 Initializing from normalized Laplacian + noise (using irlba)
## 16:27:56 Commencing optimization for 200 epochs, with 838110 positive edges
## 16:28:24 Optimization finished
DimPlot(all.anchors, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(all.anchors, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")

DimPlot(all.anchors, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "Anno1") +
DimPlot(all.anchors, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")

DimPlot(all.anchors, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "Anno2") +
DimPlot(all.anchors, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno2")

color.cnt <- scales::hue_pal()(4)[c(2,1,3,4)]
color.cnt
## [1] "#7CAE00" "#F8766D" "#00BFC4" "#C77CFF"
all.anchors$cnt <- factor(as.character(all.anchors$cnt),
                            levels = c("Stst.CTL","Stst.CKO",
                                       "CR7d.CTL","CR7d.CKO"))
DimPlot(all.anchors, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)

DimPlot(all.anchors, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1") +
DimPlot(all.anchors, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt",cols = color.cnt)

DimPlot(subset(all.anchors, subset = cnt %in% c("Stst.CTL","Stst.CKO")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[1:2]) +
DimPlot(subset(all.anchors, subset = cnt %in% c("CR7d.CTL","CR7d.CKO")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[3:4])

FeaturePlot(all.anchors, 
            reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

FeaturePlot(subset(all.anchors, subset = cnt %in% c("Stst.CTL","Stst.CKO")),
            reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

FeaturePlot(subset(all.anchors, subset = cnt %in% c("CR7d.CTL","CR7d.CKO")), 
            reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

all.anchors$sort_clusters <- factor(as.character(all.anchors$seurat_clusters),
                                    levels = c(4,6,0,2,31,14,19,13,7,8,17,
                                               5,3,9,1,18,28,11,26,
                                               10,22,21, 25,
                                               12,27,30,15, 16,24,23, 29, 20))
VlnPlot(subset(all.anchors, subset = cnt %in% c("Stst.CTL","Stst.CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "Anno1")

VlnPlot(subset(all.anchors, subset = cnt %in% c("CR7d.CTL","CR7d.CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "Anno1")

VlnPlot(subset(all.anchors, subset = cnt %in% c("Stst.CTL","Stst.CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")

VlnPlot(subset(all.anchors, subset = cnt %in% c("CR7d.CTL","CR7d.CKO")), 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")

might need to exclude mix-like C27/30 and low-quality like C31

check markers

DefaultAssay(all.anchors) <- "SCT"
DefaultAssay(all.anchors)
## [1] "SCT"
# find markers for every cluster compared to all remaining cells, with 'onlypos=F' 
Idents(all.anchors) <- "sort_clusters"

all.anchors <- PrepSCTFindMarkers(all.anchors, assay = "SCT")
options (scipen = 2)
GEX.markers.pre <- FindAllMarkers(all.anchors, only.pos = FALSE, min.pct = 0.05,
                                  assay = "SCT",
                                  test.use = "MAST",
                                  logfc.threshold = 0.25)
#GEX.markers.pre <- read.table("sn10x_integration.markers.SCT_sort.check_pre.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
## # A tibble: 256 x 7
## # Groups:   cluster [32]
##    p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene   
##    <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>  
##  1     0      1.44  0.998 0.419         0 4       Ptprt  
##  2     0      1.31  0.981 0.377         0 4       Bnc2   
##  3     0      1.21  0.994 0.514         0 4       Tshz2  
##  4     0      1.19  1     0.683         0 4       Rbfox1 
##  5     0      1.09  1     0.908         0 4       Nrg3   
##  6     0      1.06  0.998 0.792         0 4       Negr1  
##  7     0      0.944 0.938 0.483         0 4       Tox    
##  8     0      0.880 0.759 0.216         0 4       Cdc14a 
##  9     0      1.30  0.983 0.577         0 6       Galntl6
## 10     0      1.30  0.949 0.406         0 6       Grik1  
## # ... with 246 more rows

topmarkers

markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Gm",GEX.markers.pre$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                    filter(p_val_adj < 0.01) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Gm",GEX.markers.pre$gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                    filter(p_val_adj < 0.01) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
ttt = 889
ttt/60
## [1] 14.81667
ttt/64
## [1] 13.89062
ttt/65
## [1] 13.67692
pp.t120 <- list()

for(i in 1:15){
pp.t120[[i]] <- DotPlot(all.anchors, features = rev(markers.pre_t120[(60*i-59):(60*i)]), group.by = "sort_clusters")  + coord_flip() + 
           theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}

pp.t120
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

(DimPlot(subset(all.anchors,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="Stst only")) +
(DimPlot(subset(all.anchors,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="CR7d only"))

(DimPlot(subset(all.anchors,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="Stst only")) +
(DimPlot(subset(all.anchors,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="CR7d only"))

(DimPlot(subset(all.anchors,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="Stst only")) +
(DimPlot(subset(all.anchors,subset=cnt %in% c("Stst.CTL","Stst.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="Stst only"))

(DimPlot(subset(all.anchors,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "Anno1")+labs(title="CR7d only")) +
(DimPlot(subset(all.anchors,subset=cnt %in% c("CR7d.CTL","CR7d.CKO")), label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")+labs(title="CR7d only"))